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= estimated parameters for regression 
= exponent for objective function term i and variable j 
= exponent for term i, variable j, in constraint 1 
= spacecraft space debris area 

= acceleration factor of primal penalty function for constraint 1 
= spacecraft orientation factor 
= coefficient for objective function term i 
= coefficient for term i in constraint 1 
= coefficient for posyseparable term i 
= bumper material speed of sound 
= projectile diameter 

= geometric programming degree of difficulty 
= non-normalized impact velocity distribution 
= normalized impact velocity distribution 
= space debris flux 

= fraction of hyperspace for random search 
= constraint 1 
= spacecraft altitude 
= spacecraft inclination 
= number of independent variables 
= right hand side of primal constraint 1 
= wall material constant 
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m = projectile mass 

m h = number of random search points 

m = number of terms in constraint 1 

n - number of terms in objective function or number of plates (bumpers and wall) 
n . = positive integer value corresponding to variable j 

N = c umulati ve space debris flux or number of "walls" penetrated (including witness plates) 
= number of walls penetrated (normal impact) 

N t = total meteoroid flux 

p = number of constraints 

p = space debris growth rate 

P k = required confidence for random search 

p o = spacecraft probability of no penetration 
q = number of discrete variables 

r. = discrete availability factor for variable j 
s = solar flux 

S,S t = bumper/wall separation 
/ t. = bumper thickness 
t H = wall thickness 
T = mission duration 


V = projectile impact velocity 

V = maximum space debris impact velocity 
W = structure mass per unit area or weight 


a, = acceleration factor of primal penalty function for discrete constraint 1 
x = dual variable corresDonding to obiective function term i 
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«' = dual variable corresponding to term j in constraint 1 

5 ; = binary factor of primal penalty function for constraint 1 

S = first dual variable for discrete constraint of variable j 
iy 

Sj. = second dual variable for discrete constraint of variable j 

= binary factor of primal penalty function for discrete constraint 1 


e = convergence parameter for penalty function 
£j = initial exploratory step size for Hooke and Jeeves 
Ej = final exploratory step size for Hooke and Jeeves 
6 = impact angle from surface normal 

= dual objective function variable in constraint 1 
v = dual objective function 

<j> = primal penalty function 

p 1 ,p i = bumper density 


p 2 ,p„ = wall density 


= projectile mass density 
y = spacecraft inclination factor 
[ = nearest integer of quantity in brackets 

A 0 subscript denotes optimal value for a primal variable. 
A * superscript denotes optimal value for a dual variable. 
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1 INTRODUCTION 

1.1 Problem Statement 

Spacecraft designers have been concerned since the 1960’s about the effects of meteoroid 
impacts on mission safety. Recent concerns have extended to the space debris environment, 
which typically displays more massive particles than the meteoroid environment for the same 
risk level. Additionally, the higher exposure area- time product of future space missions (e.g., 
Space Station) poses a more critical design problem than current short term missions. Finally, 
the inherent uncertainties in projectile mass, velocity, density, shape, and impact angle make 
the traditional deterministic design approach impractical. 

The engineering solution to this design problem has generally been to erect a bumper or 
shield placed outboard from the spacecraft wall to disrupt/deflect the incoming projectiles. This 
passive measure has resulted in significant structural weight savings relative to a single wall 
concept with the same protective capability. The problem, then, is how to efficiently design 
these protective structures so that the bumper disrupts the projectile without posing a lethality 
problem to the wall protecting the crew and equipment. 

Spacecraft designers have a number of tools at their disposal to aid in the design process. 
These include hypervelocity impact testing, analytic impact predictors, and hydrodynamic 
codes. Perhaps the most widely accepted of these tools is impact testing, which has the advantage 
of providing actual spacecraft design verification. On the other hand, maximum test velocities 
are currently limited (8 km/sec) relative to maximum space debris (about 15 km/sec) and 
meteoroid (about 72 km/sec) velocities. Also, extensive testing is required to develop statis- 
tically significant trends for the large number of parameters associated with hypervelocity 
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impact Hydrodynamic code analysis can overcome the velocity limitation problem. However, 
this method is very computer (and time) intensive, and there is a fair amount of controversy 
involved in the selection of appropriate codes and code-specific parameters. 

Analytic impact predictors generally provide the best quick-look estimate of design 
tradeoffs. Their use is constrained by the limitations of the testing from which they are 
experimentally derived, the assumptions used in their theoretical derivation, or the regression 
analysis used in their statistical formation. However, analytic predictors may provide infor- 
mation that is clearer than that obtained from the examination of experimental results. 

The most complete way to determine the characteristics of an analytic impact predictor 
is through (nonlinear) optimization of the protective structures design problem formulated with 
the predictor of interest. Optimization techniques provide analytic or numerical solutions 
depending on the nature of the predictor, the problem formulation, and the technique used. 
1.2 Contract Purpose 

The purpose of this contract is to provide Space Station FREEDOM protective structures 
design insight through the coupling of design/material requirements, hypervelocity impact 
phenomenology, meteoroid and space debris environment sensitivities, optimization techniques 
and operations research strategies, and mission scenarios. Major findings from contract 
inception to the beginning of this study are detailed in References 100-105 and are shown below: 
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13 Study Goals 

The goals of this study are to: 


1 Develop a Monte Caiio simulation tool which will provide top level insight for Space 
SlatkaipimectivestracmieS designers* . >1 

% Develop advanced shielding concepts televant to Space Station Freedom using unique 
multiple hamper approaches* || | s l^||p! i||B | || - * 

3. Investigate projectile Shape effects on protective structuresdesign. 

The period of performance for this effort is 7-1-91 through 6-30-92. 

1.4 Study Results 

L Goal 1 was eoi^leted and i$ discussed in Section 2. A sounce listing is provided in 

AppdxA, : V P ti|§§ 

’2l?Gdid^lha*'don^bfiei4 imd is S| 
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25* Optimal Protective Structures Design is Very Sensitive toAverage Mission Altitude Above 400 km. 
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34* Shift to 15 cm Separation Results in About 33% Reduction in Protective Weight 
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2 MONTE CARLO SIMULATION DEVELOPMENT TASK 

2.1 Monte Carlo Simulation Purpose 

The purpose of this simulation is to provide a statistical tool to address and quantify 
protective structures design risks, uncertainties, and options, and to address system-level issues 
relevant to designer decision-making and possible implications. The system of initial interest 
is the structural configuration of WP01 , including the Core Module Configuration. "Grow-to" 
systems include module internal configurations and external structures (trusses, solar arrays, 
etc.) as specified in the redesign. 

Initial investigations of interest include statistical analyses of primary impacts, penetra- 
tions, and vulnerable areas. "Grow-to" investigations include interior effects, secondary ricochet 
effects, and SSF element interrelations. 

Risk considerations include environment particle velocity, impact angle, and component 
probability of impact. Uncertainty considerations include SSF IOC/FOC, particle diameter, 
mass-density, shape, and uncertainties in particle velocity and impact angle distributions. 

2.2 Monte Carlo Simulation Development Approach 

The tool development approach is to define the current SSF mission parameters and design 
configuration, and interpret the geometry mathematically using FASTGEN. The mission 
parameters drive requirements specification, including environment definitions. These con- 
siderations, combined with appropriate random number modules and the FASTGEN results, 
produce the necessary shotline time histories and intersecting body calculations. Survivability 
assessments follow and employ deterministic models for hypervelocity penetration prediction. 
Statistical assessments follow to supply answers to the questions of interest. 
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2.3 Particle Time-Arrival Process for Monte Carlo Simulation Development 

Several algorithms have been developed for the particle time-arrival process . The standard 
assumption in this area is that arrival times are Poisson distributed. This means that the inter- 
arrival times are exponentially distributed, and sorting of arrival times is not required. Mean 
data is derived from the environment flux and appropriate spacecraft areas. This algorithm 
leads to a terminating simulation defined by the mission profile. 

Realistically, however, the meteoroid and debris environments are both nonstationary 
Poisson processes, at best, since the mean arrival rates vary in time over the mission profile. 
An approximation algorithm has been developed which alters the mean arrival rate to represent 
the time period under consideration. However, this algorithm is not exact, since a period of 
high arrival rates could be neglected using a low arrival rate corresponding to the previous 
period, or vice versa. Thus, a more exact (continuous) algorithm should be developed. The 
approximating algorithm for the space debris environment is given as: 
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2. Using equations [l]-{73» develop a cumulative flux-diameter distribution: 
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In step 1, the mission start date and end date are input as reals, followed by the maximum 
and minimum particle diameters of interest and the associated particle step size. The maximum 
particle of interest is generally that particle size above which protective structures designers 
reasonably assume no liability in their considerations. For instance, there is a small but nonzero 
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probability that the Space Station will be impacted by a truck-size particle. However, designers 
could not be reasonably expected to defeat that particle size given weight constraints. Under 
current environment expectations, the choice of maximum particle diameter could reasonably 
be in the 10-1000 cm range. The choice of maximum particle diameter affects result accuracy 
and run time. The minimum particle diameter of interest is generally considered to be that size 
below which impacts are highly unlikely to lead to perforations, contamination, or any other 
measure of effectiveness which is of interest to the protective structures designer. This value 
is very much a function of the ballistic limit curves of interest Reasonable values for current 
design scenarios might be found below 0. 1 or 0.2 cm. The designer may wish to vary this value 
to investigate the sensitivity of designer measure of performance to this input. A sufficiently 
low value for this parameter would be one below which no appreciable difference resulted in 
design measure of effectiveness. It should be noted that the minimum diameter should be 
carefully selected due to its strong effect on run time and result accuracy. The diameter step 
size is also a critical factor in these two areas. This value determines the fidelity of the resulting 
flux-diameter distributions. The smaller this value is, the more accurate will be the distribution 
in a continuous sense. However, the run time increases with decreasing step size, since the flux 
distributions are being evaluated at each event time. (The requirement for re-evaluation of flux 
distributions at each event is due to the nonstationary Poisson characteristic of the space debris 
process. Generally, the diameter step size is chosen to be some fraction, say 10-20%, of the 
minimum particle diameter. The maximum, minimum, and step size values for the particle 
diameter are tuning parameters whose optimal values will vary depending on mission, con- 
figuration, and design measure of effectiveness. Thus, it is highly recommended that these 
values be adjusted before final results are displayed for each individual analysis. 
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In step 2, the flux-diameter distribution is developed at the given mission time as a function 
of the tuning parameters. This distribution is then transformed to a cumulative flux-diameter 
distribution, which is normalized to achieve a valid probability distribution from which random 
numbers may be drawn. The file size of this distribution is completely determined by the tuning 
parameters. A switch in the program input file allows the user to print out these distributions 
throughout the mission. This option should be carefully used for long missions or large dis- 
tribution sizes due to the potentially large files that would be created. 

In step 3, a uniform random variate is drawn. This variate is compared in step 4 with the 
normalized cumulative flux-diameter (probability) distribution to determine the particle 
diameter at the next event. The (absolute) value of the flux is then inverted in step 5 to determine 
the exponential (interarrival) parameter. A second uniform random variate is drawn in step 6 
for determining the interarrival time for the particle diameter determined from step 4. Step 7 
uses the exponential distribution and the results of steps 5 and 6 to determine the particle 
interairival time. In step 8, the mission clock is updated and compared with the final mission 
time in step 9. If the clock is past the end of the mission, the mission is terminated and statistics 
are gathered. Otherwise, the process reverts to step 2 for creation of new distributions and 
events. 

If independent mean and variance data for arrival rates are available, a uniform arrival 
process may be used as an alternative to Poisson arrivals. To compare this approach with the 
Poisson process, the variance may be set equal to the square of the mean. An algorithm has 
been developed for independent mean and variance data. 
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Augmentation/repair times may be modelled using a number of distributions, if this 
modelling is of interest. If mean data only for time to repair is available, an exponential service 
model may be used. If independent mean and variance data are available, the gamma, weibull, 
lognormal, or beta distributions may be appropriate. 

2.4 Simulation Status 

To date, the following items have been completed: 


A. Enveloping Geometries Established For. 

5 L Sphere: Enter Radius. v: j; ^ |fi#^ ' ' , 

X Cylinder: Enter Radius,; length 

3. Box: Enter length. Width, Height • 

B. Nonltationaiy Poisson Arrival Proeessi Algorithm Few Spade Debris ; ; i s 

H? L First Random Variate Esiablshes Point On Cumulative MuxdXameter Curve At 
Current Mission Time ♦ i 

2. Absolute Mux Is Inverted To Give Mean Intermrival Time 
• ; % Second Randon^flrtate Establishes Time Between Arrivals Using Exponential 
Distribution 

4 Cumulative Flux-Diameter Curve Is Updated For New Mission Time 
C Impact Characteristics For Space Debris 
t. Impact Velocity 

2. Impact Angle . ' . 

3, Particle DenrityAMbtis •' 

D. Look-Up Tables ... 

1. Solar Flaw (Monthly) , : 

2. Inclination Factor ; , is ' . ' . - 

3. Flux-Diameter Curves ; 

E Impact ipstory Data Including Event Time, Diameter, Density, Mass, Velocity, 
Angle ' " ' '. • 

F. Fixed Ticoe Date including Absolute Flux, Normalized Mux, Cumulative Nor- 
malizedMt4t Distributions As Functions of Diameter 

G. Geometry/Shodine integration to Include Impact Location 

H. Staridhairy Poisson AarivalProcess Algorithm- For Meteoroids ; i : rElillllliIgti 
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2.5 Selected Results 

A great many sensitivities may be generated using this simulation. Figures 2.5-1 through 

2.5- 21 provide examples of a few possible results. The following tables provide the input values 
for these example cases. (Results are for the enveloping geometry only.) Figures 2.5-22 through 

2.5- 25 show results for the current baseline configuration, shield design and mission. 

jj 1995-2004 1 1 
£ 460 km Altitude 
3. 28.5 degree Inclination 
* Configuration v | 

. 1. Spherical Enveloping Geometry 
2.: Radius — 72-5 m 
| Recorded Particles 
jjl. 0.02 -5.0 cm 

• 2. Spherical [ 

! 3. 5 % Debris Growth Rate 
*Flux Distribution Fidelity: 0.01 cm 

* ifASA TM-100471 (Debris) 
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Figure 2.5-6. Debris Impacts vs. Change in Time 
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Figure 2.5-7. Debris Impacts vs. Change in Time 
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Figure 2.5-10. Debris Impacts vs. Particle Diameter 
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Figure 2.5-11. Meteoroid Impacts vs. Particle Diameter 
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Figure 2.5-12. Debris Impacts vs. Particle Mass 
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Figure 2.5-13. Meteoroid Impacts vs. Particle Mass 
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Figure 2.5-15. Velocity Probability Distribution for 2 8.5 Degrees Inclination 
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Particle Velocity (km/sec) 

Figure 2.5-17. Meteoroid Velocity Probability Distribution 
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Figure 2.5-18. Debris Flux at Beginning and End of Mission 
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Particle Diameter (cm) 


Figure 2.5-19. Debris Flux at Beginning and End of Mission 





MULTIPLE RUN STATISTICS 


Average CPU Time = 278.61 secs ( Fidelity = 0.01 cm) 


Average CPU Time - 319.13 secs (Fidelity = 0.02 cm) 
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Figure 2.5-22. PNP VS. Mission Start Date 
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1997-2007, 28.5 deg., 

Baseline Design/Config., Combined Multibumper, 
CBD Env. With Alum. Proj. Only, Ballistic Limit. 


Figure 2.5-24. PNP VS. Average Mission Altitude 
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3 ADVANCED SHIELDING CONCEPTS TASK 

3.1 Introduction 

The development of advanced shielding concepts presented in this section includes a 
pre liminar y theoretical modification of the Wilkinson and ballistic PEN4 predictors to multiple 
bumper situations and nonlinear regression of multibumper test data from the MSFC Hyper- 
velocity Impact Test Database. 

3.2 Extension to Multiple Bumpers for Wilkinson Predictor 

A number of different approaches have been attempted to modify the Wilkinson predictor 
theoretically for multiple bumper systems. The one successful approach (physically) found 
from these approaches is given as follows: 

Modify the Wilkinson form in a product sense as: 






K 

n sf P . 

M-i ) 


0.364D 3 p,Vcos(e) D p 

for ~i — 

n p iU 

; = i 


[ 3 . 2 - 1 ] 


0.364P PpV cos(6) 


pp, 

n - 1 

n pa 

i = l 


> 1 . 


[ 3 . 2 - 2 ] 


If our goal is to minimiz e system mass per unit area subject to the total separation between first 
bumper and last wall equal to some desired value, we may write this as 


«-i 0.364D 4 p*Fcos8 


mm W - I m,+ / n _, w„, ^ 

[3.2-3] 

,=1 Ljns ? \\nmA 


V = 1 A*=» ) 


s.t. £ S,- = Sj- oT 

i = 1 

[3.2-4] 

where m i = p,^ 

[3.2-5] 
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Sr oT is the total separation between the first bumper and the wall, and n-1 is the total number of 
bumpers (n is the total number of plates). 

Under condition [3.2-2], the dual Geometric Programming objective function is given by 


* - 1 


maxv(8) = n (I/ 8;) 61 (K/8 n f 

i = 1 


WJsfe)' 


0.364D 4 ppV cos(0) 


K = 


[3.2-6] 


[3.2-7] 


X 8, = 1 

i = 1 


5, - 5, = 0, i = l,2,...,n-l 
-25„ + 8;=0, j = 1,2, ...,/i — 1 

»-l , 

m = I s ) 


[3.2-8] 

[3.2-9] 

[3.2-10] 

[3.2-11] 


Note that the degree of difficulty is 0, with 2n-2 independent variables corresponding to the n-1 
bumper areal densities and the n-1 separations. 

Equations [3.2-9] and [3.2-10] together imply 

5. = 8„ = 1/n, i = 1,2, ...,n - 1 [3.2-12] 

8; = 28 n = 2 In, j = l,2,...,n - 1 [3.2-13] 

The minimum weight and globally areal densities are given by 


2*-2 


W n = n 


O.364D 4 ppVcos(0) 


n Vn 
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0.364D 4 p*V cos(0) 
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^ $ToT ) 
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3.3 Results 

Several results using the development of Section 3.2 are given in this section. The baseline 
assumptions are a particle density of 2.8 gm/cm 3 , velocity of 9 km/sec, diameter of 1 cm, 
impacting normally into a configuration with a total bumper/wall separation of 10 cm. 

Figures 33-1 and 33-2 show how the optimal protective structures design configuration 
varies with number of bumpers for projectile diameters of 1 and 3 cm, respectively. Note that 
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for a 1 cm particle diameter, the optimal number of bumpers is 2, while for 3 cm, it is 3 bumpers. 
Also, note the significant penalty for choosing the wrong number of bumpers in these cases, as 
well as the lack of symmetry of these penalties about the optimal number of bumpers. 

Figure 33-3 shows the optimal protective structures design configuration including 
optimal number of bumpers as a function of particle diameter. Increasing particle diameter 
results in an increasing optimal number of bumpers to defeat the particle. Note the optimal 
transition regions between 1 and 2 bumpers (corresponding to particle diameters between 0.75 
and 1 cm) and 2 and 3 bumpers (corresponding to particle diameters between 2.25 and 2.5 cm). 
Also, note the very linear minimum system areal density, showing the stabilizing effect of 
increasing the number of bumpers in the configuration. 

Figure 33-4 shows the optimal protective structures design configuration including 
optimal number of bumpers as a function of particle velocity. The most striking feature of this 
trade is the relative insensitivity to velocity for a dual bumper system. 

Figure 33-5 shows the optimal protective structures design configuration as a function 
of total bumper/wall separation. As in previous studies, there is a large weight incentive for 
increasing the total separation. Furthermore, increased separation allows for more bumpers to 
disrupt the incoming particle. 

Figure 33-6 is a replica of Figure 33-5, except that the optimal individual separations 
are included. 
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Particle Density * 2.8 gm/cm3, Norma) Impact 
Particle Velocity - 9 km/sec, Diameter * 1 cm 
Total Bumper/Wall Separation * 10 cm 


Figure 3.3-1. Determining Optimal Number of Bumpers for Multibumper Wilkinson 
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Number of Bumpers 

Particle Density - 2.8 gm/cm3, Normal Impact I 

Particle Velocity - 9 Km/sec, Diameter * 3 cm I 

Total Bumper/Wall Separation- 10cm I 


Figure 3.3-2. Determining Optimal Number of Bumpers for Multibumper Wilkinson 
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Figure 3.3-3. Optimal Protective Structures Design Values vs. Particle Diameter 

(Multibumper Wilkinson) 
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Particle Density - 2.8 gm/cm3, Normal Impact 

Particle Diameter - 1 cm 

Total Bumper/Wail Separation = 10 cm 


Figure 3.3-4. Optimal Protective Structures Design Values vs. Particle Velocity 

(Multibumper Wilkinson) 
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Figure 3.3-5. Optimal Protective Structures Design Values vs. Total Bumper /Wall 

Separation (Multibumper Wilkinson) 
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Particle Density - 2.8 gm/cm3, Normal Impact 
Particle Velocity « 9 km/sec, Diameter • 1 cm 


Figure 3.3-6. Optimal Protective Structures Design Values vs. Total Bumper/Wall 

Separation (Multibumper Wilkinson) 
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3.4 Extension to Multiple Bumpers for Ballistic PEN4 Predictor 

The multiple bumper recursion equations are given by: 

7\ 

Vy- = 4100, ^ <0.4 


( T J 

V /= 498 VJ ’ F >0 ' 4 


r/ 


V 50 = 


0.67} 


( 03XlDp V' 3 

[[{—) COS(6) 


V« 3, 0C i" 2 

-2l 
Pp 


The first bumper is penetrated if 


V>V. 




The residual velocity (from the first bumper) is 

133V 2 R}p,- -(85, T^^^^cosCO) 


Vk- 


\33RpP P +R p T 1 p l /cos(Q) 


in 


The second bumper is penetrated if 




The residual velocity (from the second bumper) is 




o o ( -0.0003 125V. ^ 

1 33VgRjp p - [&S y T* 1 j/cos(0) 

133R 2 p p p +R p T 2 p 2 /cos(Q) 


in 


The third bumper is penetrated if 


[3.4-1] 

[3.4-2] 

[3.4-3] 

[3.4-4] 

[3.4-5] 

[3.4-6] 

[3.4-7] 
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3 


[3.4-8] 


The residual velocity (from the (n-l)st bumper) is 


V R = 


, , ( - 0 . 0003125V \ ,„ v nl/2 

1.33 Vjjfo, ~ [g S^J^e 

1 .33A*p, + R P T„ _ iP, _ i/ cos(0) 


[3.4-9] 


The nth bumper is penetrated if 


V * >^50 

K m - 1 *7-w 


[3.4-10] 


3.5 Advanced Shielding Task For Projectile Shatter (Multibumpers) 

The database used for regression is the MSFC Hypervelocity Impact Test Database 
Developed by the Materials & Processes Lab. Database filtering was performed by Bill Jolly 
of Sverdrup Technology Inc. Data constraints include metallic configurations with velocities 
greater than 2.5 km/sec and no MLI present. The database filtering resulted in 234 single bumper 
tests, 72 double bumper tests, and 7 triple bumper tests. A preliminary investigation using 
various posynomial regression forms was performed. The fact that Space Station Freedom has 
sufficiently low curvature in primary areas needing protection allows for the assumption of 
minimizin g system mass per unit area. The "best" intrinsically linear posynomial form resulting 
from this preliminary investigation is given by: 


n + i = jj ‘^ff s.jHfn - d'pj - 1] 
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A linear least squares analysis results in: 

K = 3.8010, a, = 1.0301, 0.3892, 03 = 0.2879, 

a 4 = 0.3288, o 5 = -0.4876, o 6 = -0.3464, 0 , = -0.2929, a, = -0.6137 [3.5-2] 

with optimal bumper scaling functions (found through search) of 

g(n - 1) = n - 1, h(n-l) = (n- l) 065 [3.5-3] 

Table 3.5-1 shows the resulting analysis of variance. 


Table 3.5-1. Analysis of Variance 


Source 

Degrees of 
Freedom 

Sum of 
Squares 

Mean Square 

F Value 

Total 

(uncorrected) 

313 

284.03 

m 

----- 

Mean 

1 

185.81 

185.81 

— 

Total 

(corrected) 

312 

9832 

03148 

----- 

Regression 

8 

2936 

3.695 

16363 

Residual 

304 

68.66 

0.2258 




The regression F value of 16363 is greater than the table value corresponding to the 
l% ^|i|f idai« e level frith 8 numerator debtees dfi&pdbift of 

:fre©d^!lTlnsital^e value iarj^ghly 231 . hypothesis that the 

slope vector of the regression plane is significantly dose to 0, and we conclude that the 

ass^piidn of a : linear relationship between the natural log vdiiesls statistical ly rigmficani. 

• ' 

^Ull — 
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The correlation coefficient is given by 

R 2 = 0.301 [3.5-4] 


R= 0.549 [3.5-5] 


The table value for the correlation coefficient ;at the 1% significancelevei, 9 total 
estimated values, and 311 degrees of freedom is well below 0,5. Since R. is greater than this 
table th^n ^Jiy pothe^ fe low ^tist^^>rre M onto this 

^Se^^ing to*a^o&cance of 0.5% and an infinite number of degrees of freedom. This 
allows ns to reject the null hypothesis that die correlation is low. 


The residual plots are shown in Figures 3.5-1 through 3.5-9. The only suspicious plot 
appears to be the one for residuals vs. number of bumpers. This correlation could be due to the 
fact that a single term posynomial does not sufficiently represent the hypervelocity impact 
phenomena over potentially complex differences between 1, 2, and 3 bumpers. 
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Particle Diameter (cm) 


5 
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Figure 3.5-1. Residuals vs Particle Diameter for Shatter Region 






Particle Density (gm/cm3) 

Figure 3.5-2. Residuals vs Particle Density for Shatter Region 
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Figure 3.5-3. Residuals vs Particle Velocity for Shatter Region 
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Cos((theta)/(n-1 )) 


Figure 3.5-4. Residuals vs Particle Impact Angle from Normal Factor for Shatter 
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Figure 3.5-7. Residuals vs Wall Mass Per Unit Area for Shatter Region 
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9.5466^'” nr^^rrF" BJ _ 

(iV + D^Cn-lf 0953 

subject to the following constraint: The total separation (first bumper to wall) is limited to a 
prespecified value 

S .t. Z S; = S ToT 

i = 1 

where m, = [3.5 — 8] 

We must determine the optimal values of the mass per unit area for the bumper(s) and wall, 
the optimal individual separations, and the minimum system mass per unit area. (S ToT is the 
total separation between the first bumper and the wall, and n-1 is the total number of bumpers 
(n is the total number of plates).) 


The problem solution is developed using Geometric Programming. The dual Geometric Pro- 
gramming objective function is given by 


max v(8) = n (1/5,) 8, (K/bS'Vh n — ]7 

i = l J s \^ToT^j 


X5, = l [3.5-11] 

; = i 


m;. 

1.1827 

‘ (n — l) a6i " ’ 

i = - 1 

[3.5-12] 

Sr 

" 1 ' 66 tls.+s;=o, 

0 n - 1) 065 ' 

j = 1,2,..., n-1 

[3.5-13] 


P 

II 

sP* 


[3.5-14] 
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These equations imply 


* 1.1827 c „ t 

5, - ^ _ jjo. 65 ^* i-l,2,...,n 1 


[3.5-15] 


o 1 

" ~ 1.1 827 (/i — 1) 0 35 + 1 


[3.5-16] 


„ 1.1827 . , „ 

o T77 i = 1,2 - 1 

* 1.1827(n-l) + (/i-l) 065 


[3.5-17] 


. 1.6647 

‘ ~ 1.1827(n — l) + (/t — l) 0 ' 65 


7-1,2,...,/* 1 


[3.5-18] 


1.6647(/i - 1) 
1.1827(n -!) + («- 1)° 


[3.5-19] 


The minimum weight and globally optimal areal densities are given by 


M827()i-1) 

1.1827(n - 1) + (? ~1)^ Y 1K7(l '~ 1) + ( ' , ~ 1)0 * i 
1.1827 J 


1.6647(n - 1) 
1.1827(w - l) + (n — l) 065 


,,K7( - 1>+( ’- ,) [ 1.1827(n-l) + (n-l) ( 
I 1.6647S ror 


0.65 

— [3.5-20] 




= [1 - (« - 


The optimal individual separations are given by 


O 101 • 1 ^ 1 

S.= J = l, 2 ,...,n - 1 

jo n -\ 
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Baseline Parameters 


The baseline parameters for an impact systemic parameter sensitivity study are given in 

Table 3.5-1. 


Table 3.5-1. Baseline Systemic Impact Parameters 


* Particle Diameter « 1 cm 
fJJ * Particle Density * 2.8 gm/cra 3 
* Particle Velocity « 5 km/sec 
f Sf Total Bumper/Wall Separation * 10 cm 
* Normal lmpact 
[£:=: * Ballistic Limit (N * 1) 


Figure 3.5-10 shows the sensitivity of the optimal protective structures design mass per 
unit area to the number of bumpers in the configuration for a particle diameter of 1 cm and 
velocity of 5 km/sec. The optimal number of bumpers is one. Figure 3.5-11 shows the same 
sensitivity except for a particle diameter of 2 cm. In this case, the optimal number of bumpers 
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is 2. Thus, the optimal number of bumpers is clearly a function of the systemic impact parameters. 

Figure 3.5-12 shows the sensitivity of optimal protective structures design variables 
(including optimal number of bumpers) to particle diameter. A transition region is found between 
1 and 1.25 cm. In this region, the optimal number of bumpers changes from 1 to 2 due to 
increases in diameter penetrability. 

Figure 3.5-13 shows the sensitivity of optimal protective structures design variables 
(including optimal number of bumpers) to particle velocity. The significant conclusion here is 
the lack of sensitivity of the design over a fairly wide velocity range (3-7.5 km/sec). 

Figure 3.5-14 shows the sensitivity of optimal protective structures design variables 
(including optimal number of bumpers) to total bumper/wall separation. A transition region is 
found between 15 and 20 cm total separation. In this region, the optimal number of bumpers 
changes from 1 to 2 due to increased spacing availability. Thus, as the allowable spacing from 
first bumper to pressure wall increases, the more incentive there is to add bumpers. 

Figure 3.5-15 shows the sensitivity of optimal protective structures design variables 
(including optimal number of bumpers) to particle impact angle from surface normal. No 
transition region is found in this sensitivity. The optimal number of bumpers remains at 1. 

Figure 3.5-16 shows the sensitivity of optimal protective structures design variables 
(including optimal number of bumpers) to particle density. A transition region is found between 
4.5 and 5 gm/cm 3 . In this region, the optimal number of bumpers changes from 1 to 2 due to 
increases in particle density penetrability. 

Figure 3.5-17 shows the sensitivity of optimal protective structures design variables 
(including optimal number of bumpers) to wall penetration factor (1 being ballistic limit). A 
transition region is found between 50% and 60% penetration. In this region, the optimal number 


of bumpers decreases from 2 to 1 due to increased allowable 
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0 1 2 3 4 5 6 

Number of Bumpers 


Particle Density - 2.8 gm/cm3, Normal Impact 
Particle Velocity - 5 km/sec, Diameter - 1 cm 
Total Bumper/Wall Separation - 10 cm, B.L. 


Figure 3.5-10. Determining Optimal Number of Bumpers for Shatter Region 
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Number of Bumpers 

Particle Density * 2.8 gm/on3, Normal Impact | 

Particle Velocity « 5 Km/sec, Diameter - 2 cm I 

Total Bumper/Wall Separation - 10 cm, B.L. I 


Figure 3.5-11. Determining Optimal Number of Bumpers for Shatter Region 
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0 


0.5 


2 


2.5 


1 1.5 

Particle Diameter (cm) 


Particle Density - 2.8 gm/cm3, Normal Impact 

Particle Velocity * 5 km/sec 

Total Bumper/Wail Separation - 10 cm, B.L 


Figure 3.5-12. Optimal Protective Structures Design Values vs. Particle Diameter 

(Shatter Region) 
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Particle Velocity (km/sec) 

Particle Density - 2.8 gm/cm3, Normal Impact a 

Particle Diameter - 1 cm 1 

Total Bumper/Wall Separation - 10 cm, B.L. I 


Figure 3.5-13. Optimal Protective Structures Design Values vs. Particle Velocity 

(Shatter Region) 
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Particle Velocity - 5 km/sec, DiametBr » 1 cm 
Total Bumper/Wall Separation - 10 cm 
Normal Impact, Ballistic Limit 


Figure 3.5-16. Optimal Protective Structures Design Values vs. Particle Density 

(Shatter Region) 
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Particle Velocity - 5 km/sec, Diameter - 1 cm 
Total Bumper/Wail Separation - 10 cm 
Normal Impact, Particle Density « 2.8 gm/cm3 


Figure 3.5-17. Optimal Protective Structures Design Values vs. Wall Penetration 

(Shatter Region) 

3.6 Multibumper Protective Structures Design Trades 

The baseline parameters and assumptions for an orbital debris analysis of multibumper 
systems is shown below. The parametric sensitivities investigated are also shown. 
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:: A* Change Board Directive Debris ]fevironment Except 
' llll Particle « Z& ss nfcxBi* A' 

II : 111 1997-2007 m^ Twsilme ■ til 1 Hill 

' > •••' °\:h " * Altitude == 398 km ■ : |j f| ' - 

* l^linauon»283deg, 

*;Total Mission PNP~&9733 (Lab;»^Q<iP^o(te witii 
§§ OipolA^ LOX Cryo,2 N* Ciycv(X9955/i^inenO -llll 

1 IbtaioMIliililiiill 1 1 

*; ' '"'V'' 1 *W r MuMbumpef | 

• "• Wilkinson/Shatter Regression (Integrated) 

% 1 -Total Bumper/WaH Separation = 10cm 
* Ballistic Limit 


\ : •" 

; * Mission Duration: 5-30 years 

* Average Mission Altitude: 200- 1000 km 

* total Mission PNP; 0.8-0.99 

* Total Debris Area; 

; : 100-1000 m 2 . > 

♦ total Bnm|»i/«rall 
j j|jjfm Separation: 5-30 cm 

Figure 3.6-1 shows the sensitivity of the optimal protective structures design mass per 
unit area to the number of bumpers in the configuration for the baseline parameters. The optimal 
number of bumpers is one. 

Figure 3.6-2 shows the sensitivity of optimal protective structures design variables (in- 
cluding optimal number of bumpers) to mission start date. A three year slip in the mission start 
date results in a 33% increase in protective structures design weight. The optimal number of 
bumpers remains constant at one for mission start dates through 2005. 
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Figure 3.6-3 shows the sensitivity of optimal protective structures design variables (in- 
cluding op timal number of bumpers) to mission duration. A transition region from one bumper 
to two is found in the 10-15 year duration range. The shape of this curve is partly reflective of 
the space debris growth rate model and partly reflective of the solar flux effect. 

Figure 3.6-4 shows the sensitivity of optimal protective structures design variables (in- 
cluding op timal number of bumpers) to average mission altitude. A transition region is found 
between 400 and 500 km altitude. In this region, the optimal number of bumpers changes from 
1 to 2 due to increased particle threat. 

Figure 3.6-5 shows the sensitivity of optimal protective structures design variables (in- 
cluding op timal number of bumpers) to total mission probability of no penetration. A transition 
from 2 bumpers to one is found in the region between 0.9733 and 0.98 PNP. This corresponds 
to element PNP’s between 0.9955 and 0.9966. 

Figure 3.6-6 shows the sensitivity of optimal protective structures design variables (in- 
cluding op timal number of bumpers) to total mission debris area. A transition region is found 
between 700 and 800 m 2 . In this region, the optimal number of bumpers changes from 1 to 2 
due to increases in particle threat size. 

Figure 3.6-7 shows the sensitivity of optimal protective structures design variables (in- 
cluding optimal number of bumpers) to total bumper/wall separation. No transition region is 
found from 5 to 30 cm total separation. An increase in total separation from 10 to 15 cm results 
in a 33% reduction in weight. 

Figures 3.6-8 through 3.6-11 show ballistic limit curves for the combined (shatter and 
0.8-Wilkinson) multibumper predictor. Seven curves are shown in each figure corresponding 
to seven different first bumper thicknesses. The baseline Space Station protective structures 
design is found behind the first bumper. The separation bet ^jy|^j^^j^j^g nd bumpers is 
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varied from 2 to 8 inches across the four curves. The transition region from projectile shatter 
to vaporization is assumed to be in the 7-8 km/sec region and is not well understood. Curves 
like these are useful in determining protective structures measures of performance for advanced 
shielding systems. 



CBD Deb. Env.(Der*s.=2.8 gm/cm3), MSHATT/MW1LK0.8 
1997-2007, 398 km, 28.5 deg., 0.9733 PNP, 603 m2 
TotaJ Bumper/Wall Separation « 10 cm, B.L. 


Figure 3.6-1. Determining Optimal Number of Bumpers 
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Mission Start Date 


CBD Deb. Env.(Dens.-2.8 gm/cm3), MSHATT/MW1LK0.8 
10 years, 398 km, 28.5 deg., 0.9733 PNP, 603 m2 
Total Bumper/Wall Separation >10 cm, B.L 


Figure 3.6-2. Optimal Protective Structures Design Values vs. Mission Start Date 
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Figure 3.6-5. Optimal Protective Structures Design Values vs. Total Mission PNP 
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SI = 5.08 cm 
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7 8 9 10 11 12 13 14 15 16 17 

Projectile Velocity (km/sec) 


tl -0.06128 cm (0.032 in.) 

— B — 

tl -0.06880 cm (0.035 In.) 
6 

tl -0.08652 cm (0.038 in.) 

— «o— - 

tl-0.10414 cm (0.041 in.) 
— * — 

tl -0.1 11 76 cm (0.044 in.) 
— — 

tl -0.1 1838 cm (0.047 In.) 

— 4 — 

tl -0.1 2700 cm (0.050 in.) 



6061-T6 Bumpers, 2219-T87 Wall. t2-0.127 cm 
twall-0.3175 cm, S2-10.16 cm, Alum. Proj. 
Normal Impact, Ballistic Limit 


Figure 3.6-9. Ballistic Limit Curve for Front Separation of 10.16 cm 
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Figure 3.6-10. Ballistic Limit Curve for Front Separation of 15.24 cm 
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Figure 3.6-11. Ballistic Limit Curve for Front Separation of 20.32 cm 


3.7 Intrinsically Nonlinear Regression For Multibumper Projectile Shatter 


The intrinsically nonlinear regression form is given by 
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To date, the addition of posynomial terms has not shown significant advantages or 
improvements over a single term posynomial for the combined multibumper database. Multiple 
term posynomial predictors, multiple databases, and separation of variables are approaches 
curiendy under investigation (but beyond the scope of this effort). 
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1* Irurinsicaily linear Fosynomial Regression C^BePerfomied to Statistically Significant 
Levels for Multiple Bumpers. " 111 ! | I 5 

2. Residual Plot vs. Number of Bumped Parameter. 

3. Other Residual Plots Seem Appropriate. ||llll|lll|ll|lll ;V - 

4. Jtasultiiqj Geometric Pmgram has ODe^^ 11 

5. Optimal Area! Densities ate Equal for Bumper(s); v > ^ Up llliil 

6. Optimal Bumpeif s) and Wall Areal Densities are Generally Not EquaL 

7. Bumper Areal Density Dormnates WaftAreal Densityfor Single Bumper Case. 

8. Wall Areal Denrily Dominates Bumper Areal Densities fcnr Multiple Bumpers. 

9. Cfctimai Individual Separations areEqual A" 

10L Optinsd Number of Bumpers Increases with Increasing Particle Diameter. (Note 
Transition Region Between d*l cm and d=T,25 cm * particle sizes.) 

11. Penalty for Selecting Wrong Number of Bumpers is Not Symmetric about Optimal 

Solutions. I | ' ' j 5 

12. Protective Structures Design Sensitivity to Velocity is Flat with Constant Optimal Number 

of Bumpers “1. ' | H , §§|| '' | , 

13. Transition from Single to DoubleBumper System is Foundfor Total Standoffs Between 

15 and 20cm. . ,s:,\ X '■ . 1 11 jf *•/ ; 

14. Impart Angle Sensitivity is Monotonically Decreasing with Constant Optimal Number 

of H ampers » 1. | f§ , - | §f ? .fs , \ 

15. Transition from Single to DoubleBumper System is Found for Particle Densities Between 

4.5 and 5 gtn/cm*. ; - 'i i 

16. Minimum SystemMass Per Unit Area is Sensitive to Wall Penetration Factor. (Transition 
From# Bumpers to I is Found Between 0.5 and 0.6.) 

17. Minimum System Mass Per Unit Area is Sensitive to the Number of Bumpers. 

IS. Three Year Schedule Slippage Results In 33% Increase in Design, 

19., Optimal Protecd^Structures Design is Very Sensitive to Mission Duration, 

20. Transition Region From l to# Bonders is Between 10 and 15 Year Durations. • . • 

2Ii| :0|rtnnal Protective Structures Design is Very: Sensitive to Average Mission Altitude 
Above 400 km. . . • , \\ . ' „ - 

22. Triutsitk» Regioi:Ffo^ 1 to 2 Bumpers is Between 400 and 500 km Altitudes. , 1 1 

23. Optimal Protective Structures Design is Very Sensitive to Mission PNP Above 0.%, 

24. Knee of die PNP Curve is Compatible With Basehne Requirement of 0.9733. 5 ; • « < 

Bumpers Is Between 0,9733 and OJ98 PNP. (0.9955 and 
Oi^966flifeanent). " . ' ! \ 


28. Optimal Protective Structures Design is Sensitive to Total BumperAVall Separation 
Between 5 and 20 cm. 

29. Knee of die Separation Curve Appears to Be Between 10 and 15 cm. : 

3©. Shift to 15 cm Separation Results in About33% Reduction in Protective Weight. 
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4 PROJECTILE SHAPE EFFECTS TASK 

SAIC developed posynomial regression techniques and combined them with posynomial 
op timization techniques for application to this area. These techniques are available for immediate 
application to the test data resulting from projectile shape effects testing. Currently, limited test 
data produces unclear results when attempts are made to correlate data from various projectile 
shapes. Results are inconclusive. Further investigation of the projectile shape effects includes 
methodologies found in sources such as "A Preliminary Investigation of Projectile Shape Effects 
In Hypervelocity Impact of a Double-Sheet Structure," by R. H. Morrison, NAS A-TN-6944, August 
1972, but will remain inconclusive until further test are performed. A summary of literature in this 
area follows. 

Burch reasoned that for multiplate systems, both projectile geometry and material play a 
significant role in terms of the characterization of hypervelocity impact damage. More specifically, 
Piekutowski provides visual confirmation of the effect that cylinder inclination at impact has on 
debris cloud features and rear wall damage severity. Shockey reports that provided the cylinder 
length exceeds both its diameter and the target plate thickness, at least part of the projectile will 
remain solid, independent of velocity. He also states that disk-type projectiles (L/D < 1) tend to be 
more completely vaporized, and that rod-like projectiles tend to remain in the solid state. Backman 
adds that pointed penetrators pierce, while blunt-shaped projectiles plug. 

Morrison considered cylinders of various L/D ratios fired at near 0 degree inclination, and 
found that cylinders were generally, and often considerably, better penetrators than spheres. (This 
is due to the lack of vaporization and shatter of the entire rod, the spike in the debris cloud, and the 
fact that the debris cloud front for cylinders travels about 14% faster than for spheres.) Specifically, 
Morrison found that on an equal mass basis, cylinders with L/D of 3/2 were the most damaging to 
double sheet structures, followed by L/D ratios of 1/2 and 1, respe^jj|i|£^MaEaf§jjycept a constant 
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L/D of 1/2, decreased L and D simultaneously, and found that a cylinder could achieve the same 
penetrability as a sphere at 1/3 the mass. Next, he kept the cylinder diameter constant, decreased 
the length systematically, and found that the same penetrability could be achieved with 1/7 the mass 
of the sphere. Counterintuitively, he discovered that for L/D ratios between 1/7 and 1/2, the lower 
ratios are more effective penetrators of 2-sheet structures, given equal mass. To quote Morrison, 
"Thus, the use of only spherical projectiles in testing double-sheet structures can be misleading. In 
fact, this practice has yielded nonconservative empirical equations now in use in the design of 
impact-resistant spacecraft structures. Although cylinders are no more typical of meteoroid shapes 
than spheres, the damage incurred by actual meteoroids is undoubtedly being underestimated by a 
significant factor." Obviously, this conclusion is even more pertinent for space debris impacts. 
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C 

c SPIES 

C Spacecraft Impact Environment Simulation 

Q ************************************************** 

C »*** Robert A. Mog ***» 

C **** SAIC Huntsville **** 

c ************************************************** 

c 

C Program to model space debris and meteoroid environments 
C for the analysis of in-orbit spacecraft survivability. 

C 

C 

c ********** REVISION JAN 14,1992 ******* 

C K FACTORS AND SURFACE AREA DIFFER POR THE BOX 

C THEREFORE ADDED SAM & SAD (SURFACE AREA METERIOD & SURFACE AREA DEBRIS) 


£***************************************** 

DOUBLE PRECISION TjDELTAT 

REAL INCL,KJ5,H,TLTF.ET ,S,SA,SAD,SAM,PVD( 150,3), 

1 DIA,DEN,MASS,VEL,ANGLE,MVEL( 1 00,2),MFL(0:1000,5) 

INTEGER N 


COMMON /ARRAYS/SOLARFLX(0: 1 50),FL(0: 1000,0:3) 

Q ***************************************************** 

C INITIALIZE PARAMETERS AND ENVELOPING GEOMETRY 

C 


CALL GET PARAMETERS(TLTFJ)MIN_D.DMAX_DJ>_D.DV, 

1 DMIN_MX>MAX_M t D_M,H,INCL, 

2 NSEEDJDISTFLAGJMFLAG) 

CALL GEOMETRY(SAD,SAMJC) 

Q ..............*..**^.**..**.***************** ******* 

C 

CALL GETSLRFLX(TI,TF) 

CALL GET_INCL_FACTOR(INCL > PSI) 

CALL DEB_VEL(INCL I PVDJ)VJMAX_PVD) 

C 

IMAX = INT((DMAX_D - DMIN_D) / D_D) 

C 

C **»* OUTPUT BEGIN AND END MISSION FLUX DISTRIBUTIONS **** 
C ***» IF IDISTFLAG IS SET **** 

C 

IF (IDISTFLAG .EQ. 1) THEN 
OPEN(UNIT=9,ST ATUS-NEW JTLE=' VELDIST.OUT) 

T = TI 

CALL FLUX_DIAMETER(SAD,DMIN_D t DMAX_D,D_D,K,HJ > SI.T.TI) 
WRITE(9,525) T 
DO 15 1=0, IMAX 

1 5 WRITE(9,530) (FL(I,II),II=0,3) 

T = TF 

CALLFLUX_DIAMETER(SAD.DMIN_D.DMAX_D,D_D.KWSI.T.TI) 
WRITE(9,525) T 
DO 16 I=0JMAX 

16 WRITE(9,530) (FL(I t n),II=0,3) 

WRITE(9,527) 

DO 17 I=1JMAX_PVD 

17 WRITE(9,520) (PVD(I,II),II=1,3) 

ENDIF 

£ ****************************************************** 

C OPEN(UNIT=8, ST ATUS-NEW 1 pTLE='SPIES_D.OUT) 

OPEN(UNIT= 1 1 ,STATUS=NEW’,FILE=TEN99.0UT) 

C 

C 

T = TI 

1 0 CALL FLUX_DLAMETER(S AD,DMIN_D,DMAX_D JC JLPSI ,T,TI) 

C 

£ *** ********** **************************************** 

C LOOK UP DIAMETER 
C 

X = RAN(NSEED) 

CALL LKUP_DIAMETER(X T DLA,IMAX,N) 

P ***************************************************** 


no noon non ~ n non non nnnnnnn nnnnn noon 
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************************** 

CALCULATE DENSITY AND MASS 

DEN = DENSITY(DIA) 

MASS = CALC_MASS<DIA) 

************************** 


*********************************************** 

LOOK UP VELOCITY AND CALCULATE IMPACT ANGLE 

X=RAN(NSEED) 

CALL LKUP.VELCPVDJMAX.PVD^VEL ANGLE) 

*********************************************** 


***************** 

REVISION JAN 14,1992 - JANEIL HILL 
SURVIVABILITY MODULE 
BASED ON VELOCITY 

*************** 

CALL SURVIVE_HIT(VEUDI AJDENANGLE) 

************************ 

CALCULATE CHANGE IN TIME 

X=RAN(NSEED) 

B = l./FL(N f l) 

DELTAT = LOG (X) * (-B) 

T = T + DELTAT 

************************ 

WRITE<8,500) T,DELTAT,DIA,DEN,MASS,VELANGLE 

IF (T.GE.TF) GOTO 150 

GOTO 10 

0 CONTINUE 
CLOSER) 

CLOSE(ll) 

******************* METEOROID ENVIRONMENT ********************* 
IF (IMFLAG .NE. 1) GOTO 200 

OPEN(UNIT=10,STATUS=’NEW’ JTLE= , SPIES_M.OUT) 

T = TI 

IMAX = INT(((DMAX_M - DMIN.M) / D_M)) 

***************************************************** 

BUILD METEOROID FLUX AND VELOCITY DISTRIBUTION TABLES 

CALL METEOROID_FLUX(SAM,H,DMIN_MJ)MAX_M 1 D_M $ MFL) 

CALL GET_METEOROID_VELOaTY (MVELJVELM AX) 

***************************************************** 


IF (IDISTFLAG EQ. 1) THEN 
WRITE(9,525) T 
DO 155 I=0JMAX 

155 WRITE(9,535) MFL(I,l).MFL(I,3),MFL(I,4),MFL<l t 5) 
WRITE(9,537) 

DO 156 1=1 JVELMAX 

156 WRITE(9,538) I,MVEL(I,1),MVEL(I,2) 

CLOSE(9) 

ENDIF 

C ********************************************* 

C LOOK UP DIAMETER 
C 

175 X=RAN(NSEED) 

CALL LKUP_MET_DIAMETER(MFLA>DIA,IMAX,N) 

Q ********************************************* 

c 

C **** DENSITY IS CONSTANT, SET MASS AND DENSITY 
MASS = MFL(N,2) 

DEN = 0.5 

IF(MFL(N,2).LE. 1E-06)DEN=2.0 

IF(MFL( N,2)XE.0.0 1 . AND.MFL<N^).GT, 1E-06)DEN= 1 .0 

Q *************************************************** 
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************************************************* 

LOOK UP VELOCITY 

X=RAN(NSEED) 

CALL LKUP_METJVEL(MVELJVELMAX,X,VEL) 

************************************************* 

**** ANGLE IS FIXED **** 

ANGLE = 90.0 
************************ 

************************************* 

CALCULATE CHANGE IN TIME 

X=RAN(NSEED) 

B = l./MFL(N t 3) 

DELTAT = LO G(X) * (-B) 

T = T + DELTAT 

Q ************************************* 

WRITEOO^IO) T»DELTAT t DIAJ)EN>lASS,VEL, ANGLE 
C 

IF (T.GE.TF) THEN 
CLOSE(10) 

CALL SORT_MERGE 
GOTO 200 
ENDDF 
C 

GOTO 175 
C 

200 CONTINUE 
STOP 

500 FORMAT(2X; D ' i 2F13.8^(2X t F6.3),2X f F10.7 t 2(2XJF6.3)) 
510 FORMAT{2X; M \2F13.8,2(2X I F6.3),2XJ ; 10.7 t 2(2XJF6.3)) 
520 FORMAT(2X,3(2X t F6.3)) 

525 FORM AT(//2X,F 13.8) 

527 FORM AT(//2X, ’Debris Velocity Distribution’) 

530 FORMAT(2X,F6.3,2X,F10.3,2(2X,F6.3)) 

535 FORMAT(2X t F6.3,2X,F10.3,3(2X,F6.3)) 

537 FORM AT(//2X /Meteoroid Velocity Distribution’) 

538 FORMAT(2X,I4,2(2X,F6.3)) 

END 


REAL FUNCTION CALC_MASS(D) 

CALC.MASS = 3.1416 * DENSITY(D) * D ** 3. / 6. 

RETURN 

END 


REAL FUNCTION DENSITY(D) 
C IF (D.LT.0.62) THEN 

C DENSITY = 4.0 

C ELSE 

C DENSITY = 2.8 / D ** 0.74 

C ENIMF 

DENSITY = Z8 

RETURN 

END 


REAL FUNCTION FLUX(KAH,PSIXS) 

REAL PHI t Fl t F2,GLG2,PSI,K 
C *** FIXED GROWTH RATE P *** 

P = 0.05 

Q=0.02 

QP=0.04 

q *************************** 

PHI = 10.**(H/200.-S/140.-1.5)/(10.**(H/200.-S/140.-1.5)+1.) 
FI = 1.22E-05 * D ** (-Z5) 

F2 = 8.1E10 * (D + 700.) ** (-6.) 

G1 = (1 + Q) ** (T-1988.) 

IF(T.GE.20 1 1 )THEN 
G1 =(1 + Q) **(23.) 

G1 = Gl*(l + QP) ** (T-2011.) 
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ENDIF 

G2 = (1 + (P *(T-1988.») 

HD=( 1 0.0* *(2.7 1 828* *((- 1 .0*(LOG(D)-0.78)* *2.0)/0.405769)))**0.5 
FLUX = K * HD*PHI * PSI * (FI * G1 + F2 * G2) 

RETURN 

END 


SUBROUTINE LKUPJXAMETER(XR>IAJMAX t N) 

COMMON / ARR AYS/SOL ARFLX(0: 1 50) RL(0: 1000,0: 3) 

I = 0 

IF (X.LT.FL(0,3)) GOTO 100 
DO 50 I = UMAX 

5 0 IF ((X.GE.FL(M,3)) .AND. (X.LT.FL(I,3))) GOTO 100 
100 CONTINUE 
DIA = FL(I,0) 

N = I 

RETURN 

END 

SUBROUTINE FLUX_DIAMETER(SAJ)MIN,DMAXJ)J<: f HJ > SI,T,TI) 

COMMON /ARRAYS/SOLARFLX(0:150)JTL(0:1000,0:3) 

REAL FLTOTJ^OB3AJDMINT>MAXT)JULPSI,T,TI 
S = SOLARFLX(INT((T-TI)*12)) 

FLTOT = 0.0 
PROB = 0.0 

IMAX = INT((DMAX-DMIN)/D) 

IF (IMAX .GT. 1000) WRITE(* *) "DATA IS TOO LARGE FOR \ 

I DEBRIS FLUX ARRAY: FL’ 

DO 100 I = 0.IMAX 

DIAMETER = DMIN + I*D 
FL(I,0) = DIAMETER 

FL(U) = SA * FLUX(KJ>IAMETER,HRSI t T,S) 

FLTOT = FLTOT + FL(U) 

100 CONTINUE 

DO 110 I = 0,IMAX 
FL(U) = FL(I,1) / FLTOT 
PROB = PROB + FL(U) 

FL(I,3) = PROB 
110 CONTINUE 

C *** INDICES FOR FL ARE AS FOLLOWS: 

C 0 - DIAMETER / 1 - IMPACTS/YEAR / 2 - NORMALIZED / 3 - CUMULATIVE 

Q ***************************************************************** 

RETURN 

END 

SUBROUTINE GET_INCL_FACTOR(INCLJPSI) 

REAL INCL 

(X*EN(UNTT= 12JTLE= , FLUXFAC.DAT f STAT T US= OLD’ READONLY) 

DC = 0 

DO 10 WHILE (DC XT. INCL) 

1 0 READ(12 t *RND=20) IX RSI 

GOTO 30 

2 0 WRITER*) ERROR READING FLUXFAC.DAT 

3 0 CLOSE(12) 

RETURN 

END 

SUBROUTINE GETSLRFLX(TI,TF) 

COMMON / ARR AYS/SOLAR FLX(0 : 1 50) RL(0: 1000, 0:3) 

OPEN(UNTT= 1 2RILE='SOLARl .FLX . ST ATUS^’OLD’ READONLY) 

II = INT((TI- 1988.)* 1 2.) 

DO 10 1=1,11 

10 READ(12,*) 

II = INT((TF-TI)*12.) 

DO 20 1=0 M 
READ(12 *RRR=25) S 
SOLARFLX(I) = S 
20 CONTINUE 
GOTO 30 

25 WRITE( V) 1 ERROR READING SOLAR 1 .FLX 
3 0 CLOSE(12) 

RETURN 

END 
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SUBROUTINE GEOMETRY (S AD,S AMJC) 

REAL R^ t W,H,PI,K 
PI = 3.141593 
5 WRITE(MOO) 

READ(M10) I 
IF(I.LT.l .OR. I.GT.3) THEN 
WRITE(* *) 1 Invalid Choice. ’ 

GOTO 5 
ENDIF 

GOTO (10,20,30), I 
1 0 WRITE(\210) 

READ(*,215) R 
SAD = 4*PI*R*R 
SAM = 4*PI*R*R 
K = Z7 
GOTO 50 

2 0 WRITE(*,220) 

READ(*,215) R 
WRITE(*,225) 

READ(*,215) L 

SAD = 2* PI * R * (L + R) 

SAM = 2* PI * R * (L + R) 

C 

C *** NEED TO CALCULATE K*** 

C 

GOTO 50 

30 WRITE(*,230) 

READ(*,215) L 
WRITE(*,232) 

READ(*^15) W 
WRITE(* t 235) 

READ<*,215) H 
K1 = 2.6 
K2= 1.6 

SAD= K1*W*H + 2*K2*L*H 
SAM = 2*L*W + 2*W*H + 2*L*H 
K= 1. 

C 

C *** WITH BOX A DOUBLE K FACTOR IS NEEDED *** 

C FOR NOW IT IS HARD WIRED IN K1 = 2.6 (FRONT) K2 = 1.6 (SIDE) 

C 

50 WRITE(*,300) SAD.SAM 

RETURN 

100 FORMAT(2X, ’Input type of enveloping geometry:’/! Ox ,'l - Sphere’/ 
1 10X/2 - CylinderyiOX/3 - Box.’/2X,’Choice?\$) 

110 FORMAT(I2) 

210 FORM AT(4X,' Sphere :76X,’Enter Radius (m) :’,$) 

215 FORMAT(F7.2) 

220 FORM AT(4X,’Cy Under :'/6X, ’Enter Radius (m) :\$) 

225 FORM AT(6X, ’Enter Length (m) :’,$) 

230 FORM AT(4X, ’Box :76X/Enter Length (m) :’,$) 

232 FORM AT(6X ,’Enter Width (m) :’,$) 

235 FORMAT(6X,Enter Height (m) :’£) 

300 FORMAT(4X, Surface Area is ’,F19.2,' m2,F19.2,’ m2’) 

END 

SUBROUTINE DEB_VEL(XINCLXPVJ)VJMAX) 

REAL XPV(150,3) 

YG=250.0 

YF=0.0 

YC=0.0125 

YE=0.55+0.005*(XINCL-30.0) 

YH=1.0-0.0000757*(XINCL-60.0)**2.0 

YA=2.5 

YB=0.3 

YD= 1.3-0.01 *(XINCL- 30.0) 

YV0=7.7 

IF(XINCL .LE. 60.0)THEN 
YB=0.5 
YG=18.7 

YVO=7.25+0.015*(XINCL-30.0) 

ENDIF 

EF(XINCL L£. 80.0 .AND. XINCL .GT. 60.0)THEN 
YB=0.5-0.0 1 *(XINCL-60.0) 

YG=18.7+0.0289*(XINCL-60.0)**3.0 
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ENDIF 

IF(XINCL .GT. 100.0)THEN 

YC=0.01 25+0.00 125*(XINCL- 1 00.0) 

ENDIF 

IF(XINCL .LE. 50.0)THEN 

YF=0.3+0.0008*(XINCL-50.0)**2.0 

ENDIF 

IF(XINCL .GT. 50.0 .AND. XINCL .LE. 80.0)THEN 
YF=0.3-0.01 *(XINCL-50.0) 

ENDIF 

XSUMIV=0.0 

IVMAX=1 

rv=i 

xv=i. 

100 XPV(IV,1)=XV 

XPV(TV,2)=Y G*2.7 1 83** (- 1.0*((X V-Y A* YV0)/(YB* YVO))* *2.0) 
XPV(IV^)=:XPV(IV,2)+YF*2.7183**(- 1 .0*((XV-YD*YV0)/ 

1 (YE*YV0))**2.O) 

XPV(TV,2)=XPV(IV,2)*(2.0*XV*YV0-XV**2.0) 

XPV(IV;2)=XPV(IV t 2)+YH*YC*(4.0*XV*YV0-XV**2.0) 

IF(XPV(IV,2) .LE. 0.000)THEN 
XPV(IV.2)=0.000 
IVMAX=IV 
GO TO 150 
ENDIF 

XSUMI V =XSUMI V +XPV(I V,2) 

IV=IV+1 
XV=XV+DV 
GOTO 100 
150 PROB = 0. 

DO 200 I=1*IVMAX 

XPV(L2)=XPV(I.2)/XSUMIV 
PROB = PROB + XPV(U) 

XPV(I,3) = PROB 
200 CONTINUE 

q ************** INDICES ARE: **************** 

C 1 - VELOCITY 2 - NORMALIZED 3 - CUMULATIVE 

Q ******************************************** 

IMAX = IVMAX 
C 

RETURN 

END 

SUBROUTINE LKUP_VEL(XPV t IMAXJLXV f A) 

REAL XPV(150,3) 

XV =XPV(1,1) 

DO 10 I = UMAX-1 

IF ((R.GE.XPV(I,3)) .AND. (R.LT.XPV(I+13))) THEN 
XV = XPV(I+1,1) 

GOTO 20 
ENDIF 

10 CONTINUE 
2 0 CONTINUE 

A = ACOSD(XV/l 5.4) 

RETURN 

END 

SUBROUTINE GET_PARAMETERS(TI,TF t DMIN_D,DM AX_D,D_D,DV ,DMIN_M, 
1 DMAX_M,D_M,H t INCL,NSEEDJDISTFLAG,IMFLAG) 

REAL TLTF,DMIN JDMAX,D,DV.H,INCL 
INTEGER NSEED 

OPEN(UNIT= 1 2JFILE='SPIES99. IN ,STATUS=OLD’ .READONLY) 

READ(12,*) TI 
READ(12,*) TF 
READ(12,*) DMIN_D 
READ(12,*) DMAX_D 
READ(12 *) D_D 
READ(12,*) DV 
READ(12 t *) DMIN.M 
READ(12,*) DMAX.M 
READ(12.*) D_M 
READ(12.*) H 
READ(12*) INCL 
READ( 12.*) NSEED 
READ(12,*) IDISTFLAG 
READ(12,*) IMFLAG 
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CLOSE(12) 

RETURN 

END 

SUBROUTINE METEOROID_FLUX(SA,HJ)MIN,DMAXJ),MFL) 

REAL DMIN,DMAX^,DLAMETERJLOGM,MASS,H,G,NET A.MFLTOT, 

1 MFL(0: 1000,5) 

C *** CALCULATION FOR G USES 6378 KM RADIUS FOR EARTH *** 

G » 1.0+6478./(6378.+H) 

NETA = (1 + COS(ASIN(6478. / (6378. + H))» / 2, 

IMAX = INT((DMAX-DMIN)/D) 

IF (IMAX .GT. 1000) WRITER ,♦) ’MUST RE-DIMENSION METEOROID 
1 FLUX ARRAY. DATA IS TOO LARGE.’ 

MFLTOT = 0. 

PROB = 0. 

DO 10 1=0, IMAX 

C *** MASS CALCULATION USES AVERAGE DENSITY OF .5 G/ CM3 *** 

DIAMETER = DMIN + I*D 
DEN = 0.5 

IF(MFL(N^)LE.1E-06)DEN=ZO 

IF(MFL(N,2).LE.0.01.AND.MFL(N,2).GT.1E-06)DEN=1.0 
MASS = 3.141593 * DEN * DIAMETER ** 3.0 / 6.0 
MFL(U) = DIAMETER 
MFL(I,2) = MASS 
OX LOGM = LOG10(MASS) 

OX IF (MASS LE. .000001) THEN 

err MFL(I,3) = 10**(-14.339 - 1.584*LOGM - 0.063 * LOGM ** 2) 

OX ELSE 

OX MFL(I,3) = 10**(- 14.37 - 1.213 * LOGM) 

OX ENDIF 

C0=3.156E07 

C1=2.2E03 

C2=15. 

C3=1.3E-09 

C4=1.0E11 

C5=1.0E27 

C6=1.3E-16 

C7=1.0E06 

MFL(I,3)=C6/(MASS+C7*MASS**2.0)**0.85 

MFL(I,3)=MFL(I,3)+C3/(MASS+C4*MASS**Z0+C5*MASS**4.0)**0.36 

MFL(I,3)=MFL(I,3)+1.0/(C1*MASS**0.306+C2)**4.38 

MFL(I,3)=C0*MFL(I,3) 

C *** MULTIPLY BY G.NETA, AND NO. OF SEC/YEAR *** 

MFL(I,3) = SA * MFL(I,3) * G * NETA * 31536000. 

MFLTOT = MFLTOT + MFL(I,3) 

1 0 CONTINUE 

DO 20 I=0,IMAX 

MFL(I,4) = MFL(I,3) / MFLTOT 
PROB = PROB + MFL(I,4) 

MFL(1 ,5) = PROB 
20 CONTINUE 

C *** INDICES ARE : ********************** 

C 1 - DIAMETER 2 - MASS 3 - IMPACTS/YEAR 
C 4 - NORMALIZED 5 - CUMULATIVE 
RETURN 
END 

SUBROUTINE GET METEOROID JVELOCTTY (MVELJMAX) 

REAL MVEL( 100,2 ),PVEL,PVTOT .PVCUM 
OPEN(UNIT=8JTLE= , METVEL.IN’,STATUS= , OLD\READONLY) 

PVCUM = 0. 

PVTOT = 0. 

1 0 READ(8 ♦,END=100)I,PVEL 
GX MVEL(I,1) = PVEL 

IF(I.LT.l 1)MVEL(I,1)=0.0 
IF(I.GE. 1 1 . AND.I.LE. 1 6)M VEL(I,1 )=0. 1 1 2 
IF(I.GT.16.AND.I.LT.55)MVEL(I,1)=3.328E05/(FLOAT(I)**5.34) 
IF(I.GE.55.AND. I.LE.72)M VEL(I, 1 )= 1 .695E- 04 
PVTOT = PVTOT + MVEL(I,1) 

GOTO 10 
100 IMAX = I 

DO 110 1=1 .IMAX 

MVELG.l) = MVELG.l) / PVTOT 
PVCUM = PVCUM + MVEL(I,1) 

MVEL(U) = PVCUM 
1 10 CONTINUE 
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CLOSE(8) 

RETURN 

BSD 

SUBROUTINE LXUP_MET_VEL(MVEL,IMAXP,VEL) 

REAL VEL.MVEL( 100,2) 

VEL = 1.0 

DO 10 I = l.IMAX-1 

IF ((R.GE.MVEUU)) .AND. (R.LT.MVEL(I+U))) THEN 
VEL = FLOAT© 

GOTO 20 
ENDIF 

10 CONTINUE 
20 CONTINUE 
RETURN 
END 

SUBROUTINE LKUP_MET_DLAMETER(MFUX,DIA JMAX.N) 
REAL MFL(0: 1000.5) 

I = 0 

IF (X.LT.MFL(0.5)) GOTO 100 
DO 50 I = UMAX 

50 IF «X.GEMFL(I-14)) AND. (XiT.MFL(I,5))) GOTO 100 
100 CONTINUE 

DIA = MFL(I.l) 

N = I 

RETURN 

END 

SUBROUTINE SORT MERGE 

CHARACTER*58 DDAT A2,M DAT A2.DD AT A 1 * 3 ,MD AT A 1 * 3 
REALTD.TM 

OPENCUNTT=8^TATUS= , OLD'pILE=’SPIES_D.OUTJ?EADONLY) 
(X’EN(UNTr=9STATUS='OLD'pILE='SPIES_M.OUT READONLY) 
OPEN(UNIT=10,STATUS= , NEW’pILE='SPIES_BOTH.OLT) 
READ(8,500pND= 100) DDATA1.TD.DDATA2 
READ(9,500JEND=200) MDATA1.TM.MDATA2 
50 DO 60 WHILE (TD PE. TM) 

WRITE( 10,500) DDATA1.TD.DDATA2 
READ(8,500pND= 100) DDAT A 1 .TD.DD AT A2 
60 CONTINUE 

DO 70 WHILE (TD .GT. TM) 

WRITE( 10,500) MDAT A 1 .TM.MDATA2 
READ(9,500pND=200) MDATA1.TM.MDATA2 
7 0 CONTINUE 
GOTO 50 

100 READ(9,500,END=300) MDATA1.TM.MDATA2 
WRITE<10,500) MDATA1.TM.MDATA2 
GOTO 100 

200 READ(8,500,END=300) DDATA1 , TD.DD ATA2 
WRITE(10,500) DDATA1.TD.DDATA2 
GOTO 200 
300 CLOSE(8) 

CLOSE(9) 

CLOSE(10) 

RETURN 

500 FORMAT(2X.A3,F13.8,A58) 

END 


SUBROUTINE SURVTVE_HIT( VEL, DIA.DEN, THETA) 

REAL VEL 

£******************** 

C SURVIVABILITY MODULE JAN 14,1992 
C JANEIL HILL 
C VEL VELOCITY OF PARTICLE 
C BIGN PERFERATION LEVEL 
C DIA DIAMETER OF PARTICLE 
C EEN PARTICLE DENSITY 

C THETA ANGLE OF IMPACT 
C SI ARRAY OF SEPERATION OF BUMPER 

C RHOI ARRAY OF DENSITY OF BUMPER 

C THI BUMPER THICKNESS 

C N NUMBER OF BUMPERS 

Q************ ********** 

C FOR NOW ALL THE VARIABLES THAT WILL BE READ IN ARE HARDWIRED HERE 
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REAL SI(10).RHOI(10).THI(10)LN 

LN = .401 
N = 2 

THI(l) = .127 
THI(2) = .3175 
SI(1) = 10.16 
RHOI(l) = 2.71 
RHO!(2) = 2.81 


IF (VEL X.E 7.5) THEN 
TEMPI = 1.0 
TEMP2 = 1.0 
DO I=1,(N-1) 

TEMPI = TEMPI * SI(I) 

TEMP2 = TEMP2 *(RHOI(I) * THI(I)) 

END DO 

BIGN =( ( 3.801*DIA** 1 .0301 * DEN*».3892 * VEL**.2879 * 

& (COSD(THETAAN-l))**.3288)) / 

& (TEMPI **(.4826 / (N-l)**.65) * TEMP2**(.3464/(N-1)**.65) * 

& (RHOI(N)*THI(N))**.2929 * (N-l)**.6137) ) - 1. 

WRITE(11,*) BIGN 
ELSE 

C **** IF VELOCITY > 7.5) 

TEMPI = 1.0 
TEMP2 = 1.0 
DO I=1,(N-1) 

TEMPI = TEMPI * SI(I)**2. 

TEMP2 = TEMP2 *(RHOI(I) * THI(D) 

END DO 

RMPUA = (1.25*DIA * DEN) / TEMP2 
IF (RMPUA GT. 1.0000) THEN 

TN = (.364 * (1.25*DIA)**4. *DEN**2. * VEL * COSD(THET A)) / 
& (LN * TEMPI * TEMP2 * RHOI(N)) 

ELSE 

C **** RMPUA <= 1 

TN = (.364 *(1.25*DIA)**3.» DEN* VEL * COSD(THETA) ) / 

& (LN * TEMPI * RHOI(N) ) 

END IF 

IF (TN .GE. THI(N)) THEN 
WRITE(11 *) ' 1.500’ 

ELSE 

WRITE (11,*) ’ 0.000’ 

ENDIF 

ENDIF 

RETURN 

END 



